#ifndef dijet_H
#define dijet_H
#include "utility.h"
#include <stdio.h>
#include <iostream>
#include "mt2_bisect.h"

using namespace std;
using namespace mt2_bisect;


//****************************************************************************
extern "C" pair<long double, long double> M1OverM2_mass_pair(int a, int b,
		int c, int d, vector<vector<int> > matrixtoadd, int kp,
		vector<TLorentzVector> vec_part) {
	//****************************************************************************

	TLorentzVector sum1, sum2;
	long double deltares, mavg, M1, M2;
	pair<long double, long double> res;

	//	cout<< matrixtoadd[0][0]<< " "<< matrixtoadd[1][0]<< " " << matrixtoadd[2][0]<<" "<< matrixtoadd[3][0]<<endl;

	sum1 = vec_part[matrixtoadd[a][kp]] + vec_part[matrixtoadd[b][kp]];

	sum2 = vec_part[matrixtoadd[c][kp]] + vec_part[matrixtoadd[d][kp]];

	M1 = sqrt(fabs(sum1 * sum1));
	M2 = sqrt(fabs(sum2 * sum2));

//	deltares = fabs(1 - M2 / M1);
	deltares = fabs((M1 - M2) / (M1 + M2));

	mavg = (M1 + M2) / 2.;

	res.first = deltares;
	res.second = mavg;

	//cout<< res.first << " "<< res.second << endl;

	return res;

}

//****************************************************************************
extern "C" pair<long double, long double> M1OverM2_DeltaCMS_pair(int a, int b,
		int c, int d, vector<vector<int> > matrixtoadd, int kp,
		vector<TLorentzVector> vec_part) {
	//****************************************************************************

	TLorentzVector sum1, sum2;
	long double deltares, mavg, M1, M2, DeltaCMS, pt1, pt2;
	pair<long double, long double> res;

	//	cout<< matrixtoadd[0][0]<< " "<< matrixtoadd[1][0]<< " " << matrixtoadd[2][0]<<" "<< matrixtoadd[3][0]<<endl;

	sum1 = vec_part[matrixtoadd[a][kp]] + vec_part[matrixtoadd[b][kp]];

	sum2 = vec_part[matrixtoadd[c][kp]] + vec_part[matrixtoadd[d][kp]];

	M1 = sqrt(fabs(sum1 * sum1));
	M2 = sqrt(fabs(sum2 * sum2));

	deltares = fabs(1 - M2 / M1);
	mavg = (M1 + M2) / 2.;

	pt1 = vec_part[matrixtoadd[a][kp]].Pt();
	pt2 = vec_part[matrixtoadd[b][kp]].Pt();

//	pt1 = vec_part[matrixtoadd[c][kp]].Pt();
//	pt2 = vec_part[matrixtoadd[d][kp]].Pt();

	DeltaCMS = pt1 + pt2 - mavg;

	res.first = deltares;
	res.second = DeltaCMS;

	//cout<< res.first << " "<< res.second << endl;

	return res;

}

//****************************************************************************
long double delta_rel_res_of_4jets(vector<int> v,
		vector<vector<int> > matrixtoadd, int kp,
		vector<TLorentzVector> vec_part) {
	//****************************************************************************

	TLorentzVector sum1, sum2;
	long double deltares, mavg, M1, M2;
	long double res;

	//	cout<< matrixtoadd[0][0]<< " "<< matrixtoadd[1][0]<< " " << matrixtoadd[2][0]<<" "<< matrixtoadd[3][0]<<endl;

	int a, b, c, d;

	a = v[0];
	b = v[1];
	c = v[2];
	d = v[3];

	sum1 = vec_part[matrixtoadd[a][kp]] + vec_part[matrixtoadd[b][kp]];

	sum2 = vec_part[matrixtoadd[c][kp]] + vec_part[matrixtoadd[d][kp]];

	M1 = sqrt(fabs(sum1 * sum1));
	M2 = sqrt(fabs(sum2 * sum2));

//	deltares = fabs(1 - M2 / M1);
	deltares = fabs((M1 - M2) / (M1 + M2));

	res = deltares;

	//cout<< res.first << " "<< res.second << endl;

	return res;

}

//****************************************************************************
extern "C" pair<long double, long double> delta_mass_pair(int a, int b, int c,
		int d, vector<vector<int> > matrixtoadd, int kp,
		vector<TLorentzVector> vec_part) {
	//****************************************************************************

	TLorentzVector sum1, sum2;
	long double deltares, mavg, M1, M2;
	pair<long double, long double> res;

	//	cout<< matrixtoadd[0][0]<< " "<< matrixtoadd[1][0]<< " " << matrixtoadd[2][0]<<" "<< matrixtoadd[3][0]<<endl;

	sum1 = vec_part[matrixtoadd[a][kp]] + vec_part[matrixtoadd[b][kp]];

	sum2 = vec_part[matrixtoadd[c][kp]] + vec_part[matrixtoadd[d][kp]];

	M1 = sqrt(fabs(sum1 * sum1));
	M2 = sqrt(fabs(sum2 * sum2));

	deltares = fabs(M1 - M2);
	mavg = (M1 + M2) / 2.;

	res.first = deltares;
	res.second = mavg;

	//cout<< res.first << " "<< res.second << endl;

	return res;

}

//****************************************************************************
extern "C" pair<long double, pair<long double, long double> > delta_dy_y_res_pair_pair(
		int a, int b, int c, int d, vector<vector<int> > matrixtoadd, int kp,
		vector<TLorentzVector> vec_part) {
	//****************************************************************************

	TLorentzVector sum1, sum2;
	long double deltares, M1, M2, dyres, y2, y1;
	pair<long double, pair<long double, long double> > res;
	pair<long double, long double> deltaYy1;

	//	cout<< matrixtoadd[0][0]<< " "<< matrixtoadd[1][0]<< " " << matrixtoadd[2][0]<<" "<< matrixtoadd[3][0]<<endl;

	sum1 = vec_part[matrixtoadd[a][kp]] + vec_part[matrixtoadd[b][kp]];

	sum2 = vec_part[matrixtoadd[c][kp]] + vec_part[matrixtoadd[d][kp]];

	M1 = sqrt(fabs(sum1 * sum1));
	M2 = sqrt(fabs(sum2 * sum2));

//	deltares = fabs(1 - M2 / M1);
	deltares = fabs((M1 - M2) / (M1 + M2));

	//cout<< res.first << " "<< res.second << endl;

	y1 = sum1.Rapidity();
	y2 = sum2.Rapidity();
	dyres = fabs(y1 - y2);

	deltaYy1.first = dyres;
	deltaYy1.second = y1;

	res.first = deltares;
	res.second = deltaYy1;

	//cout<< res.first << " "<< res.second << endl;

	return res;

}

//****************************************************************************
extern "C" pair<long double, long double> delta_b2b_pair(int a, int b, int c,
		int d, vector<vector<int> > matrixtoadd, int kp,
		vector<TLorentzVector> vec_part) {
	//****************************************************************************

	TLorentzVector sum1, sum2, sum;
	long double deltares, M1, M2, dyres, y2, y1;
	pair<long double, long double> res;

	//	cout<< matrixtoadd[0][0]<< " "<< matrixtoadd[1][0]<< " " << matrixtoadd[2][0]<<" "<< matrixtoadd[3][0]<<endl;

	TVector3 boost;

	sum = vec_part[matrixtoadd[a][kp]] + vec_part[matrixtoadd[b][kp]]
			+ vec_part[matrixtoadd[c][kp]] + vec_part[matrixtoadd[d][kp]];

	boost = sum.BoostVector();
	//		cout << Ycm << endl;
	//	Ycm = 0.5 * log((sum.E() + sum.Pz()) / (sum.E() - sum.Pz()));

	//		cout << "~~~~~~~~~~" << endl;
	//		cout << sum.Vect().Mag() << endl;
	//		//			cout << (vec_part[vectorstoadd[0]] + vec_part[vectorstoadd[1]]).Vect().Mag() << endl;
	//			cout << BoostZ(sum, Ycm).Vect().Mag() << endl;
	TLorentzVector temp;
	temp = sum;
	temp.Boost(-boost);

	//cout << " Total 3-momentum of the 4b in the CoM frame " << temp.Px() << " "
	//	<< temp.Py() << " " << temp.Pz() << " " << endl;

	sum1 = vec_part[matrixtoadd[a][kp]] + vec_part[matrixtoadd[b][kp]];
	sum2 = vec_part[matrixtoadd[c][kp]] + vec_part[matrixtoadd[d][kp]];

	M1 = sqrt(fabs(sum1 * sum1));
	M2 = sqrt(fabs(sum2 * sum2));

//	deltares = fabs(1 - M2 / M1);
	deltares = fabs((M1 - M2) / (M1 + M2));

//	cout << "before the boost the masses are " << M1 << " " << M2 << endl;

	// ******** boost back in the CoM ************
	sum1.Boost(-boost);
	sum2.Boost(-boost);

	M1 = sqrt(fabs(sum1 * sum1));
	M2 = sqrt(fabs(sum2 * sum2));

	deltares = fabs(1 - M2 / M1);
//	cout << "after the boost the masses are " << M1 << " " << M2 << endl;
	//cout<< res.first << " "<< res.second << endl;

	y1 = sum1.Rapidity();
	y2 = sum2.Rapidity();

	res.first = deltares;
	res.second = fabs(2. * (y1 + y2) / (y1 - y2));

	// if y1=-y2 as for 2 resonances in the CoM this should be zero.

	//cout<< res.first << " "<< res.second << endl;

	return res;

}

//****************************************************************************
extern "C" pair<long double, long double> delta_dy_res_pair(int a, int b, int c,
		int d, vector<vector<int> > matrixtoadd, int kp,
		vector<TLorentzVector> vec_part) {
	//****************************************************************************

	TLorentzVector sum1, sum2;
	long double deltares, M1, M2, dyres, y2, y1;
	pair<long double, long double> res;

	//	cout<< matrixtoadd[0][0]<< " "<< matrixtoadd[1][0]<< " " << matrixtoadd[2][0]<<" "<< matrixtoadd[3][0]<<endl;

	sum1 = vec_part[matrixtoadd[a][kp]] + vec_part[matrixtoadd[b][kp]];

	sum2 = vec_part[matrixtoadd[c][kp]] + vec_part[matrixtoadd[d][kp]];

	M1 = sqrt(fabs(sum1 * sum1));
	M2 = sqrt(fabs(sum2 * sum2));

//	deltares = fabs(1 - M2 / M1);
	deltares = fabs((M1 - M2) / (M1 + M2));

	//cout<< res.first << " "<< res.second << endl;

	y1 = sum1.Rapidity();
	y2 = sum2.Rapidity();
	dyres = fabs(y1 - y2);

	res.first = deltares;
	res.second = dyres;

	//cout<< res.first << " "<< res.second << endl;

	return res;

}

////****************************************************************************
//extern "C" long double deltaR_res_4into2(int a, int b, int c, int d,
//		vector<vector<int> > matrixtoadd, int kp,
//		vector<TLorentzVector> vec_part) {
//	//****************************************************************************
//
//	TLorentzVector pa, pb, pc, pd;
//	long double deltares, dR1, dR2, detares, res;
//
//	//	cout<< matrixtoadd[0][0]<< " "<< matrixtoadd[1][0]<< " " << matrixtoadd[2][0]<<" "<< matrixtoadd[3][0]<<endl;
//
//	pa = vec_part[matrixtoadd[a][kp]];
//	pb = vec_part[matrixtoadd[b][kp]];
//	pc = vec_part[matrixtoadd[c][kp]];
//	pd = vec_part[matrixtoadd[d][kp]];
//
//	dR1 = pa.Vect().DeltaR(pb.Vect());
//	dR2 = pc.Vect().DeltaR(pd.Vect());
//
//	deltares = fabs(dR1 - 1) + fabs(dR2 - 1);
//
//	res = deltares;
//
//	//cout<< res.first << " "<< res.second << endl;
//
//	return res;
//
//}


//****************************************************************************
extern "C" pair<long double, long double> delta_dCMS_res_pair(int a, int b,
		int c, int d, vector<vector<int> > matrixtoadd, int kp,
		vector<TLorentzVector> vec_part) {
	//****************************************************************************

	TLorentzVector sum1, sum2;
	long double deltares, M1, M2, mavg, detares, eta2, eta1, pta, ptb, ptc, ptd,
			dCMS1, dCMS2, dCMSavg, dCMSmin;
	pair<long double, long double> res;

	//	cout<< matrixtoadd[0][0]<< " "<< matrixtoadd[1][0]<< " " << matrixtoadd[2][0]<<" "<< matrixtoadd[3][0]<<endl;

	pta = vec_part[matrixtoadd[a][kp]].Pt();
	ptb = vec_part[matrixtoadd[b][kp]].Pt();
	ptc = vec_part[matrixtoadd[c][kp]].Pt();
	ptd = vec_part[matrixtoadd[d][kp]].Pt();

	sum1 = vec_part[matrixtoadd[a][kp]] + vec_part[matrixtoadd[b][kp]];

	sum2 = vec_part[matrixtoadd[c][kp]] + vec_part[matrixtoadd[d][kp]];

	M1 = sqrt(fabs(sum1 * sum1));
	M2 = sqrt(fabs(sum2 * sum2));

	mavg = (M1 + M2) / 2.;

//	deltares = fabs(1 - M2 / M1);
	deltares = fabs((M1 - M2) / (M1 + M2));

	//cout<< res.first << " "<< res.second << endl;

	eta1 = sum1.Eta();
	eta2 = sum2.Eta();
	detares = fabs(eta1 - eta2);

	dCMS1 = pta + ptb - mavg;
	dCMS2 = ptc + ptd - mavg;

	dCMSavg = (dCMS1 + dCMS2) / 2.;
	dCMSmin = min(dCMS1, dCMS2);

	res.first = deltares;
	res.second = dCMSmin;

	//cout<< res.first << " "<< res.second << endl;

	return res;

}

//****************************************************************************
extern "C" pair<long double, long double> delta_deta_res_pair(int a, int b,
		int c, int d, vector<vector<int> > matrixtoadd, int kp,
		vector<TLorentzVector> vec_part) {
	//****************************************************************************

	TLorentzVector sum1, sum2;
	long double deltares, M1, M2, detares, eta2, eta1;
	pair<long double, long double> res;

	//	cout<< matrixtoadd[0][0]<< " "<< matrixtoadd[1][0]<< " " << matrixtoadd[2][0]<<" "<< matrixtoadd[3][0]<<endl;

	sum1 = vec_part[matrixtoadd[a][kp]] + vec_part[matrixtoadd[b][kp]];

	sum2 = vec_part[matrixtoadd[c][kp]] + vec_part[matrixtoadd[d][kp]];

	M1 = sqrt(fabs(sum1 * sum1));
	M2 = sqrt(fabs(sum2 * sum2));

//	deltares = fabs(1 - M2 / M1);
	deltares = fabs((M1 - M2) / (M1 + M2));

	//cout<< res.first << " "<< res.second << endl;

	eta1 = sum1.Eta();
	eta2 = sum2.Eta();
	detares = fabs(eta1 - eta2);

	res.first = deltares;
	res.second = detares;

	//cout<< res.first << " "<< res.second << endl;

	return res;

}

//****************************************************************************
extern "C" pair<long double, long double> delta_dphi_res_pair(int a, int b,
		int c, int d, vector<vector<int> > matrixtoadd, int kp,
		vector<TLorentzVector> vec_part) {
	//****************************************************************************

	TLorentzVector sum1, sum2;
	long double deltares, M1, M2, dphires;
	pair<long double, long double> res;

	//	cout<< matrixtoadd[0][0]<< " "<< matrixtoadd[1][0]<< " " << matrixtoadd[2][0]<<" "<< matrixtoadd[3][0]<<endl;

	sum1 = vec_part[matrixtoadd[a][kp]] + vec_part[matrixtoadd[b][kp]];

	sum2 = vec_part[matrixtoadd[c][kp]] + vec_part[matrixtoadd[d][kp]];

	M1 = sqrt(fabs(sum1 * sum1));
	M2 = sqrt(fabs(sum2 * sum2));

//	deltares = fabs(1 - M2 / M1);
	deltares = fabs((M1 - M2) / (M1 + M2));

	//cout<< res.first << " "<< res.second << endl;

	dphires = sum1.Vect().DeltaPhi(sum2.Vect()); // get deltaphi between v1 and v2

//	cout << "***************" << endl;
//	cout << sum1.Phi() << " " << sum2.Phi() << endl;
//	cout << "delta_phi("<< a << b << c << d << ")" << dphires << endl;
//	cout << "***************" << endl;

	res.first = deltares;
	res.second = dphires;

	//cout<< res.first << " "<< res.second << endl;

	return res;

}

//****************************************************************************
extern "C" pair<long double, long double> delta_deta_intraresonance_pair(int a,
		int b, int c, int d, vector<vector<int> > matrixtoadd, int kp,
		vector<TLorentzVector> vec_part) {
	//****************************************************************************

	TLorentzVector sum1, sum2;
	long double deltares, dEtaAvg, M1, M2, deltaEta1, deltaEta2, etaA, etaB,
			etaC, etaD;
	pair<long double, long double> res;

	//	cout<< matrixtoadd[0][0]<< " "<< matrixtoadd[1][0]<< " " << matrixtoadd[2][0]<<" "<< matrixtoadd[3][0]<<endl;

	sum1 = vec_part[matrixtoadd[a][kp]] + vec_part[matrixtoadd[b][kp]];

	sum2 = vec_part[matrixtoadd[c][kp]] + vec_part[matrixtoadd[d][kp]];

	M1 = sqrt(fabs(sum1 * sum1));
	M2 = sqrt(fabs(sum2 * sum2));

	etaA = vec_part[matrixtoadd[a][kp]].Eta();
	etaB = vec_part[matrixtoadd[b][kp]].Eta();
	etaC = vec_part[matrixtoadd[c][kp]].Eta();
	etaD = vec_part[matrixtoadd[d][kp]].Eta();

	//	deltares = fabs(M1 - M2);
//	deltares = fabs(1 - M2 / M1);
	deltares = fabs((M1 - M2) / (M1 + M2));

	dEtaAvg = (fabs(etaA - etaB) + fabs(etaC - etaD)) / 2.;

	res.first = deltares;
	res.second = dEtaAvg;

	//cout<< res.first << " "<< res.second << endl;

	return res;

}

//****************************************************************************
extern "C" pair<long double, long double> delta_dR_intraresonance_pair(int a,
		int b, int c, int d, vector<vector<int> > matrixtoadd, int kp,
		vector<TLorentzVector> vec_part) {
	//****************************************************************************

	TLorentzVector sum1, sum2;
	long double deltares, dRavg, M1, M2, dR1, dR2;
	pair<long double, long double> res;

	//	cout<< matrixtoadd[0][0]<< " "<< matrixtoadd[1][0]<< " " << matrixtoadd[2][0]<<" "<< matrixtoadd[3][0]<<endl;

	sum1 = vec_part[matrixtoadd[a][kp]] + vec_part[matrixtoadd[b][kp]];

	sum2 = vec_part[matrixtoadd[c][kp]] + vec_part[matrixtoadd[d][kp]];

	M1 = sqrt(fabs(sum1 * sum1));
	M2 = sqrt(fabs(sum2 * sum2));

	dR1 = vec_part[matrixtoadd[a][kp]].Vect().DeltaR(
			vec_part[matrixtoadd[b][kp]].Vect()); // get angle between v1 and v2

	dR2 = vec_part[matrixtoadd[c][kp]].Vect().DeltaR(
			vec_part[matrixtoadd[d][kp]].Vect()); // get angle between v1 and v2

	//	deltares = fabs(M1 - M2);
//	deltares = fabs(1 - M2 / M1);
	deltares = fabs((M1 - M2) / (M1 + M2));

	dRavg = (dR1 + dR2) / 2.;

	res.first = deltares;
	res.second = dRavg;

	//cout<< res.first << " "<< res.second << endl;

	return res;

}

//****************************************************************************
extern "C" pair<long double, long double> delta_dphi_intraresonance_pair(int a,
		int b, int c, int d, vector<vector<int> > matrixtoadd, int kp,
		vector<TLorentzVector> vec_part) {
	//****************************************************************************

	TLorentzVector sum1, sum2;
	long double deltares, dphiavg, M1, M2, dphi1, dphi2;
	pair<long double, long double> res;

	//	cout<< matrixtoadd[0][0]<< " "<< matrixtoadd[1][0]<< " " << matrixtoadd[2][0]<<" "<< matrixtoadd[3][0]<<endl;

	sum1 = vec_part[matrixtoadd[a][kp]] + vec_part[matrixtoadd[b][kp]];

	sum2 = vec_part[matrixtoadd[c][kp]] + vec_part[matrixtoadd[d][kp]];

	M1 = sqrt(fabs(sum1 * sum1));
	M2 = sqrt(fabs(sum2 * sum2));

	dphi1 = vec_part[matrixtoadd[a][kp]].Vect().DeltaPhi(
			vec_part[matrixtoadd[b][kp]].Vect()); // get angle between v1 and v2

	dphi2 = vec_part[matrixtoadd[c][kp]].Vect().DeltaPhi(
			vec_part[matrixtoadd[d][kp]].Vect()); // get angle between v1 and v2

	//	deltares = fabs(M1 - M2);
//	deltares = fabs(1 - M2 / M1);
	deltares = fabs((M1 - M2) / (M1 + M2));

	dphiavg = ( fabs(dphi1) + fabs(dphi2) ) / 2.;

	res.first = deltares;
	res.second = dphiavg;

	//cout<< res.first << " "<< res.second << endl;

	return res;

}

//****************************************************************************
extern "C" long double b2bness_from_y_resonances_bestpair_boost(
		vector<vector<int> > matrixtoadd, vector<TLorentzVector> vec_part,
		vector<TRootLHEFParticle*> part) {
	//****************************************************************************
	// compute the pai r of float (mass difference , average deltaphi of the  dijet)
	// and returns the average deltaphi of the pairing with smaller delta massa.
	//printmatrix(matrixtoadd);

	//TLorentzVector sum4b;

	long double b2b, deltamin, res;
	vector < pair<long double, long double> > delta_b2b;
	b2b = -100;

	if (matrixtoadd.size() == 4) {

		// delta_dphi_res_pair uses the relative best pairing
		delta_b2b.push_back(
				delta_b2b_pair(0, 1, 2, 3, matrixtoadd, 0, vec_part));
		delta_b2b.push_back(
				delta_b2b_pair(0, 2, 1, 3, matrixtoadd, 0, vec_part));
		delta_b2b.push_back(
				delta_b2b_pair(0, 3, 2, 1, matrixtoadd, 0, vec_part));

		deltamin = delta_b2b[0].first;
		b2b = delta_b2b[0].second;

		for (int k = 0; k < delta_b2b.size(); k++) {
			//			cout << delta_dy_y[k].first << " "<< delta_dy_y[k].second << endl;

			if (1. / delta_b2b[k].first > 1 / deltamin) {
				b2b = delta_b2b[k].second;
			}

		}

		res = b2b;

		//cout << res << endl;
		return res;

	} else {

		cout
				<< " the mass of the best pairing can be computed only on 4 objects"
				<< endl;
		exit(1);

	}
}

//****************************************************************************
extern "C" long double b2bness_from_y_resonances_bestpair_no_boost(
		vector<vector<int> > matrixtoadd, vector<TLorentzVector> vec_part,
		vector<TRootLHEFParticle*> part) {
	//****************************************************************************
	// compute the pai r of float (mass difference , average deltaphi of the  dijet)
	// and returns the average deltaphi of the pairing with smaller delta massa.
	//printmatrix(matrixtoadd);

	TLorentzVector sum4b;

	long double delta_y_bestpair, deltamin, y4b, y1, res;
	vector < pair<long double, pair<long double, long double> > > delta_dy_y;
	delta_y_bestpair = -100.;
	y1 = 1000.;

	if (matrixtoadd.size() == 4) {

		// delta_dphi_res_pair uses the relative best pairing
		delta_dy_y.push_back(
				delta_dy_y_res_pair_pair(0, 1, 2, 3, matrixtoadd, 0, vec_part));
		delta_dy_y.push_back(
				delta_dy_y_res_pair_pair(0, 2, 1, 3, matrixtoadd, 0, vec_part));
		delta_dy_y.push_back(
				delta_dy_y_res_pair_pair(0, 3, 2, 1, matrixtoadd, 0, vec_part));

		deltamin = delta_dy_y[0].first;
		delta_y_bestpair = (delta_dy_y[0].second).first;
		y1 = (delta_dy_y[0].second).second;

		for (int k = 0; k < delta_dy_y.size(); k++) {
			//			cout << delta_dy_y[k].first << " "<< delta_dy_y[k].second << endl;

			if (1. / delta_dy_y[k].first > 1 / deltamin) {
				delta_y_bestpair = (delta_dy_y[k].second).first;
				y1 = (delta_dy_y[k].second).second;
			}

		}

		sum4b = vec_part[matrixtoadd[0][0]] + vec_part[matrixtoadd[1][0]]
				+ vec_part[matrixtoadd[2][0]] + vec_part[matrixtoadd[3][0]];

		y4b = sum4b.Rapidity();

		res = min(fabs(delta_y_bestpair - fabs(2. * (y1 - y4b))),
				fabs(delta_y_bestpair - fabs(2. * (y1 + y4b))));

		//cout << res << endl;
		return res;

	} else {

		cout
				<< " the mass of the best pairing can be computed only on 4 objects"
				<< endl;
		exit(1);

	}
}

//****************************************************************************
extern "C" long double delta_rapidity_resonances_bestpair(
		vector<vector<int> > matrixtoadd, vector<TLorentzVector> vec_part,
		vector<TRootLHEFParticle*> part) {
	//****************************************************************************
	// compute the pai r of float (mass difference , average deltaphi of the  dijet)
	// and returns the average deltaphi of the pairing with smaller delta massa.
	//printmatrix(matrixtoadd);

	long double delta_y_bestpair, deltamin;
	vector < pair<long double, long double> > delta_massa;
	delta_y_bestpair = -100.;

	if (matrixtoadd.size() == 4) {

		// delta_dphi_res_pair uses the relative best pairing
		delta_massa.push_back(
				delta_dy_res_pair(0, 1, 2, 3, matrixtoadd, 0, vec_part));
		delta_massa.push_back(
				delta_dy_res_pair(0, 2, 1, 3, matrixtoadd, 0, vec_part));
		delta_massa.push_back(
				delta_dy_res_pair(0, 3, 2, 1, matrixtoadd, 0, vec_part));

		deltamin = delta_massa[0].first;
		delta_y_bestpair = delta_massa[0].second;

		for (int k = 0; k < delta_massa.size(); k++) {
			//			cout << delta_massa[k].first << " "<< delta_massa[k].second << endl;

			if (1. / delta_massa[k].first > 1 / deltamin) {
				delta_y_bestpair = delta_massa[k].second;
			}

		}

		// cout << Mmin << endl;
		return delta_y_bestpair;

	} else {

		cout
				<< " the mass of the best pairing can be computed only on 4 objects"
				<< endl;
		exit(1);

	}
}

//****************************************************************************
extern "C" long double delta_eta_resonances_bestpair(
		vector<vector<int> > matrixtoadd, vector<TLorentzVector> vec_part,
		vector<TRootLHEFParticle*> part) {
	//****************************************************************************
	// compute the pai r of float (mass difference , average deltaphi of the  dijet)
	// and returns the average deltaphi of the pairing with smaller delta massa.
	//printmatrix(matrixtoadd);

	long double delta_eta_bestpair, deltamin;
	vector < pair<long double, long double> > delta_massa;
	delta_eta_bestpair = -100.;

	if (matrixtoadd.size() == 4) {

		// delta_dphi_res_pair uses the relative best pairing
		delta_massa.push_back(
				delta_deta_res_pair(0, 1, 2, 3, matrixtoadd, 0, vec_part));
		delta_massa.push_back(
				delta_deta_res_pair(0, 2, 1, 3, matrixtoadd, 0, vec_part));
		delta_massa.push_back(
				delta_deta_res_pair(0, 3, 2, 1, matrixtoadd, 0, vec_part));

		deltamin = delta_massa[0].first;
		delta_eta_bestpair = delta_massa[0].second;

		for (int k = 0; k < delta_massa.size(); k++) {
			//			cout << delta_massa[k].first << " "<< delta_massa[k].second << endl;

			if (1. / delta_massa[k].first > 1 / deltamin) {
				delta_eta_bestpair = delta_massa[k].second;
			}

		}

		// cout << Mmin << endl;
		return delta_eta_bestpair;

	} else {

		cout
				<< " the mass of the best pairing can be computed only on 4 objects"
				<< endl;
		exit(1);

	}
}

//****************************************************************************
extern "C" long double delta_CMS_bestpair_ratio(
		vector<vector<int> > matrixtoadd, vector<TLorentzVector> vec_part,
		vector<TRootLHEFParticle*> part) {
	//****************************************************************************
	// compute the pai r of float (mass difference , average deltaphi of the  dijet)
	// and returns the average deltaphi of the pairing with smaller delta massa.
	//printmatrix(matrixtoadd);

	long double deltaCMS_bestpair, deltamin;
	vector < pair<long double, long double> > delta_massa;
	deltaCMS_bestpair = -100.;

	if (matrixtoadd.size() == 4) {

		// delta_dphi_res_pair uses the relative best pairing
		delta_massa.push_back(
				delta_dCMS_res_pair(0, 1, 2, 3, matrixtoadd, 0, vec_part));
		delta_massa.push_back(
				delta_dCMS_res_pair(0, 2, 1, 3, matrixtoadd, 0, vec_part));
		delta_massa.push_back(
				delta_dCMS_res_pair(0, 3, 2, 1, matrixtoadd, 0, vec_part));

		deltamin = delta_massa[0].first;
		deltaCMS_bestpair = delta_massa[0].second;

		for (int k = 0; k < delta_massa.size(); k++) {
			//			cout << delta_massa[k].first << " "<< delta_massa[k].second << endl;

			if (1. / delta_massa[k].first > 1 / deltamin) {
				deltaCMS_bestpair = delta_massa[k].second;
			}

		}

		// cout << Mmin << endl;
		return deltaCMS_bestpair;

	} else {

		cout
				<< " the mass of the best pairing can be computed only on 4 objects"
				<< endl;
		exit(1);

	}
}

//****************************************************************************
extern "C" long double delta_phi_resonances_bestpair_ratio(
		vector<vector<int> > matrixtoadd, vector<TLorentzVector> vec_part,
		vector<TRootLHEFParticle*> part) {
	//****************************************************************************
	// compute the pai r of float (mass difference , average deltaphi of the  dijet)
	// and returns the average deltaphi of the pairing with smaller delta massa.
	//printmatrix(matrixtoadd);

	long double delta_phi_bestpair, deltamin;
	vector < pair<long double, long double> > delta_massa;
	delta_phi_bestpair = -100.;

	if (matrixtoadd.size() == 4) {

		// delta_dphi_res_pair uses the relative best pairing
		delta_massa.push_back(
				delta_dphi_res_pair(0, 1, 2, 3, matrixtoadd, 0, vec_part));
		delta_massa.push_back(
				delta_dphi_res_pair(0, 2, 1, 3, matrixtoadd, 0, vec_part));
		delta_massa.push_back(
				delta_dphi_res_pair(0, 3, 2, 1, matrixtoadd, 0, vec_part));

		deltamin = delta_massa[0].first;
		delta_phi_bestpair = delta_massa[0].second;

		for (int k = 0; k < delta_massa.size(); k++) {
			//			cout << delta_massa[k].first << " "<< delta_massa[k].second << endl;

			if (1. / delta_massa[k].first > 1 / deltamin) {
				delta_phi_bestpair = delta_massa[k].second;
			}

		}

		// cout << Mmin << endl;
		return delta_phi_bestpair;

	} else {

		cout
				<< " the mass of the best pairing can be computed only on 4 objects"
				<< endl;
		exit(1);

	}
}


//****************************************************************************
extern "C" long double delta_eta_intradijet_bestpair(
		vector<vector<int> > matrixtoadd, vector<TLorentzVector> vec_part,
		vector<TRootLHEFParticle*> part) {
	//****************************************************************************
	// compute the pai r of float (mass difference , average deltaphi of the  dijet)
	// and returns the average deltaphi of the pairing with smaller delta massa.
	//printmatrix(matrixtoadd);

	long double delta_eta_bestpair, deltamin;
	vector < pair<long double, long double> > delta_massa;
	delta_eta_bestpair = -100.;

	if (matrixtoadd.size() == 4) {

		delta_massa.push_back(
				delta_deta_intraresonance_pair(0, 1, 2, 3, matrixtoadd, 0,
						vec_part));
		delta_massa.push_back(
				delta_deta_intraresonance_pair(0, 2, 1, 3, matrixtoadd, 0,
						vec_part));
		delta_massa.push_back(
				delta_deta_intraresonance_pair(0, 3, 2, 1, matrixtoadd, 0,
						vec_part));

		deltamin = delta_massa[0].first;
		delta_eta_bestpair = delta_massa[0].second;

		for (int k = 0; k < delta_massa.size(); k++) {
			//			cout << delta_massa[k].first << " "<< delta_massa[k].second << endl;

			if (1. / delta_massa[k].first > 1 / deltamin) {
				delta_eta_bestpair = delta_massa[k].second;
			}

		}

		// cout << Mmin << endl;
		return delta_eta_bestpair;

	} else {

		cout
				<< " the mass of the best pairing can be computed only on 4 objects"
				<< endl;
		exit(1);

	}
}

//****************************************************************************
extern "C" long double delta_R_intradijet_bestpair(
		vector<vector<int> > matrixtoadd, vector<TLorentzVector> vec_part,
		vector<TRootLHEFParticle*> part) {
	//****************************************************************************
	// compute the pai r of float (mass difference , average deltaphi of the  dijet)
	// and returns the average deltaphi of the pairing with smaller delta massa.
	//printmatrix(matrixtoadd);

	long double delta_R_bestpair, deltamin;
	vector < pair<long double, long double> > delta_massa;
	delta_R_bestpair = -100.;

	if (matrixtoadd.size() == 4) {

		delta_massa.push_back(
				delta_dR_intraresonance_pair(0, 1, 2, 3, matrixtoadd, 0,
						vec_part));
		delta_massa.push_back(
				delta_dR_intraresonance_pair(0, 2, 1, 3, matrixtoadd, 0,
						vec_part));
		delta_massa.push_back(
				delta_dR_intraresonance_pair(0, 3, 2, 1, matrixtoadd, 0,
						vec_part));

		deltamin = delta_massa[0].first;
		delta_R_bestpair = fabs(delta_massa[0].second);

		for (int k = 0; k < delta_massa.size(); k++) {
			//			cout << delta_massa[k].first << " "<< delta_massa[k].second << endl;

			if (1. / delta_massa[k].first > 1 / deltamin) {
				delta_R_bestpair = delta_massa[k].second;
			}

		}

		// cout << Mmin << endl;
		return delta_R_bestpair;

	} else {

		cout
				<< " the mass of the best pairing can be computed only on 4 objects"
				<< endl;
		exit(1);

	}
}

//****************************************************************************
extern "C" long double delta_phi_intradijet_bestpair(
		vector<vector<int> > matrixtoadd, vector<TLorentzVector> vec_part,
		vector<TRootLHEFParticle*> part) {
	//****************************************************************************
	// compute the pai r of float (mass difference , average deltaphi of the  dijet)
	// and returns the average deltaphi of the pairing with smaller delta massa.
	//printmatrix(matrixtoadd);

	long double delta_phi_bestpair, deltamin;
	vector < pair<long double, long double> > delta_massa;
	delta_phi_bestpair = -100.;

	if (matrixtoadd.size() == 4) {

		delta_massa.push_back(
				delta_dphi_intraresonance_pair(0, 1, 2, 3, matrixtoadd, 0,
						vec_part));
		delta_massa.push_back(
				delta_dphi_intraresonance_pair(0, 2, 1, 3, matrixtoadd, 0,
						vec_part));
		delta_massa.push_back(
				delta_dphi_intraresonance_pair(0, 3, 2, 1, matrixtoadd, 0,
						vec_part));

		deltamin = delta_massa[0].first;
		delta_phi_bestpair = fabs(delta_massa[0].second);

		for (int k = 0; k < delta_massa.size(); k++) {
			//			cout << delta_massa[k].first << " "<< delta_massa[k].second << endl;

			if (1. / delta_massa[k].first > 1 / deltamin) {
				delta_phi_bestpair = fabs(delta_massa[k].second);
			}

		}

		// cout << Mmin << endl;
		return delta_phi_bestpair;

	} else {

		cout
				<< " the mass of the best pairing can be computed only on 4 objects"
				<< endl;
		exit(1);

	}
}

//****************************************************************************
extern "C" long double deltaM_bestpair(vector<vector<int> > matrixtoadd,
		vector<TLorentzVector> vec_part, vector<TRootLHEFParticle*> part) {
	//****************************************************************************
	// compute the total invariant mass of N particles
	// particle involved are those in the event with IDs listed in  the vector vectorstoadd
	//printmatrix(matrixtoadd);

	long double Mmin, deltamin;
	vector < pair<long double, long double> > delta_massa;
	Mmin = -100.;

	if (matrixtoadd.size() == 4) {

		delta_massa.push_back(
				delta_mass_pair(0, 1, 2, 3, matrixtoadd, 0, vec_part));
		delta_massa.push_back(
				delta_mass_pair(0, 2, 1, 3, matrixtoadd, 0, vec_part));
		delta_massa.push_back(
				delta_mass_pair(0, 3, 2, 1, matrixtoadd, 0, vec_part));

		deltamin = delta_massa[0].first;
		Mmin = delta_massa[0].second;

		for (int k = 0; k < delta_massa.size(); k++) {
			//cout << delta_massa[k].first << endl;

			if (1. / delta_massa[k].first > 1 / deltamin) {
				deltamin = delta_massa[k].first;
				Mmin = delta_massa[k].second;
			}

		}

		// cout << Mmin << endl;
		return deltamin;

	} else {

		cout
				<< " the mass of the best pairing can be computed only on 4 objects"
				<< endl;
		exit(1);

	}
}

//****************************************************************************
extern "C" long double minv_bestpair(vector<vector<int> > matrixtoadd,
		vector<TLorentzVector> vec_part, vector<TRootLHEFParticle*> part) {
	//****************************************************************************
	// compute the total invariant mass of N particles
	// particle involved are those in the event with IDs listed in  the vector vectorstoadd
	//printmatrix(matrixtoadd);

	long double Mmin, deltamin;
	vector < pair<long double, long double> > delta_massa;
	Mmin = -100.;

	if (matrixtoadd.size() == 4) {

		delta_massa.push_back(
				delta_mass_pair(0, 1, 2, 3, matrixtoadd, 0, vec_part));
		delta_massa.push_back(
				delta_mass_pair(0, 2, 1, 3, matrixtoadd, 0, vec_part));
		delta_massa.push_back(
				delta_mass_pair(0, 3, 2, 1, matrixtoadd, 0, vec_part));

		deltamin = delta_massa[0].first;
		Mmin = delta_massa[0].second;

		for (int k = 0; k < delta_massa.size(); k++) {
			//cout << delta_massa[k].first << " "<< delta_massa[k].second << endl;

			if (1. / delta_massa[k].first > 1 / deltamin) {
				deltamin = delta_massa[k].first;
				Mmin = delta_massa[k].second;
			}

		}

		//cout << Mmin << endl;
		return Mmin;

	} else {

		cout
				<< " the mass of the best pairing can be computed only on 4 objects"
				<< endl;
		exit(1);

	}
}

//****************************************************************************
extern "C" long double ratio_bestpair_ratio(vector<vector<int> > matrixtoadd,
		vector<TLorentzVector> vec_part, vector<TRootLHEFParticle*> part) {
	//****************************************************************************
	// compute the total invariant mass of N particles
	// particle involved are those in the event with IDs listed in  the vector vectorstoadd
	//printmatrix(matrixtoadd);

	long double Mmin, deltamin;
	vector < pair<long double, long double> > delta_massa;
	Mmin = -100.;

	if (matrixtoadd.size() == 4) {

		delta_massa.push_back(
				M1OverM2_mass_pair(0, 1, 2, 3, matrixtoadd, 0, vec_part));
		delta_massa.push_back(
				M1OverM2_mass_pair(0, 2, 1, 3, matrixtoadd, 0, vec_part));
		delta_massa.push_back(
				M1OverM2_mass_pair(0, 3, 2, 1, matrixtoadd, 0, vec_part));

		deltamin = delta_massa[0].first;
		Mmin = delta_massa[0].second;

		for (int k = 0; k < delta_massa.size(); k++) {
			//cout << delta_massa[k].first << " "<< delta_massa[k].second << endl;

			if (1. / delta_massa[k].first > 1 / deltamin) {
				deltamin = delta_massa[k].first;
				Mmin = delta_massa[k].second;
			}

		}

		//cout << Mmin << endl;
		return deltamin;

	} else {

		cout
				<< " the mass of the best pairing can be computed only on 4 objects"
				<< endl;
		exit(1);

	}
}

//****************************************************************************
extern "C" long double minv_bestpair_ratio(vector<vector<int> > matrixtoadd,
		vector<TLorentzVector> vec_part, vector<TRootLHEFParticle*> part) {
	//****************************************************************************
	// compute the total invariant mass of N particles
	// particle involved are those in the event with IDs listed in  the vector vectorstoadd
	//printmatrix(matrixtoadd);

	long double Mmin, deltamin;
	vector < pair<long double, long double> > delta_massa;
	Mmin = -100.;

	if (matrixtoadd.size() == 4) {

		delta_massa.push_back(
				M1OverM2_mass_pair(0, 1, 2, 3, matrixtoadd, 0, vec_part));
		delta_massa.push_back(
				M1OverM2_mass_pair(0, 2, 1, 3, matrixtoadd, 0, vec_part));
		delta_massa.push_back(
				M1OverM2_mass_pair(0, 3, 2, 1, matrixtoadd, 0, vec_part));

		deltamin = delta_massa[0].first;
		Mmin = delta_massa[0].second;

		for (int k = 0; k < delta_massa.size(); k++) {
			//cout << delta_massa[k].first << " "<< delta_massa[k].second << endl;

			if (1. / delta_massa[k].first > 1 / deltamin) {
				deltamin = delta_massa[k].first;
				Mmin = delta_massa[k].second;
			}

		}

		//cout << Mmin << endl;
		return Mmin;

	} else {

		cout
				<< " the mass of the best pairing can be computed only on 4 objects"
				<< endl;
		exit(1);

	}
}

//****************************************************************************
extern "C" long double howmany_pairings_better_than_75pc(
		vector<vector<int> > matrixtoadd, vector<TLorentzVector> vec_part,
		vector<TRootLHEFParticle*> part) {
	//****************************************************************************
	// compute the total invariant mass of N particles
	// particle involved are those in the event with IDs listed in  the vector vectorstoadd
	//printmatrix(matrixtoadd);

	vector < pair<long double, long double> > delta_massa;
	int howmany;
	howmany = 0;
	long double ratiomax;
	ratiomax = 0.075;

	if (matrixtoadd.size() == 4) {

		delta_massa.push_back(
				M1OverM2_mass_pair(0, 1, 2, 3, matrixtoadd, 0, vec_part));
		delta_massa.push_back(
				M1OverM2_mass_pair(0, 2, 1, 3, matrixtoadd, 0, vec_part));
		delta_massa.push_back(
				M1OverM2_mass_pair(0, 3, 2, 1, matrixtoadd, 0, vec_part));

		for (int k = 0; k < delta_massa.size(); k++) {
			//cout << delta_massa[k].first << " "<< delta_massa[k].second << endl;

			if (delta_massa[k].first < ratiomax) {
				howmany++;
			}

		}

		//cout << Mmin << endl;
		return howmany;

	} else {

		cout
				<< " the mass of the best pairing can be computed only on 4 objects"
				<< endl;
		exit(1);

	}
}

#endif
